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ABSTRACT 

We discuss the behavior of gas dynamic flows which are perturbations of a 
uniform stream in terms of information transfer across artificial 
(computational) boundaries remote from the source of disturbance. A set of 
boundary conditions are derived involving vorticity, entropy, and pressure- 
velocity relationships derived from bicharacteristic equations. 
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1. INTRODUCTION 


A recurring frustration in Computational Fluid Dynamics is the apparent 
difficulty of giving numerical expression to very simple statements. Typical 
of this situation is the specification, in aerodynamic problems, that the flow 
is uniform at large distances. The problem is caused by the fact that the 
outer limit of the computational domain never is truly at infinity, because it 
can always be reached by a numerical signal. Therefore, merely specifying 
uniform conditions on an outer boundary results in an overconstrained problem; 
signals which do reach the outer boundary are liable to be reflected from it 
and may completely corrupt the interior solution. 

There has been a long search for effective absorbing boundary conditions, 
but none so far has found universal acceptance, and many practical codes make 
use of empirical procedures. In this note, we indicate the fallacy in three 
current practices and advocate a new procedure which may be less objectionable 
and can be applied to unsteady flow in any number of dimensions. 

We remark that there is no problem when the flow is supersonic at 
infinity. It is then both simple and correct to prescribe everything at 
inflow and nothing at outflow. Upper and lower boundaries can be treated as 
rigid walls remote enough that reflected waves do not impinge on the region of 
interest. Our treatment, therefore, concentrates on the subsonic case where 
the difficulties are twofold. The boundary conditions must lead to a well- 
posed problem, so as to avoid the instabilities associated with overcon- 
straint, and they should also be an accurate statement of the physics so that 
they can be applied at fairly small distances. In the absence of rigorous 
analysis (which is difficult, see [1] for a recent review), we hope that the 
second property will imply the first. We derive physically correct equations 
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which express the passage of different kinds of information across the 
boundary. If these are differenced so that the numerical information used is 
taken from the proper domain of dependence, we assume that stability will be 
assured. 


2. SOME UNSOUND PRACTICES 

2.1 SPECIFY FLOW DIRECTION AT INFLOW 

Consider Figure 1. The proposal is to set v = 0 on AB. This means 

that AB cannot contribute to the circulation u»ds around ABCD. The 
circulation should be independent of the shape or distance of the integration 
contour (provided the flow is inviscid and all shocks are inside the 
contour). It is understandable, therefore, that with this method lifting 
forces are usually underestimated [2]. A cure, applicable to steady two- 
dimensional aerofoil flows, is to match the direction to that found in a far- 
field analytic solution whose circulation would produce the lift measured on 
the aerofoil at each iteration. No such procedure is available for three- 
dimensional or unsteady flows. 


2.2 SPECIFY PRESSURE AT OUTFLOW 

This is an essentially empirical procedure which often works quite 
well. However, it would not be applicable in a three-dimensional flow with 
shed vorticity. In steady incompressible flow, for example, we should satisfy 
Bernoulli's equation 
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1 2 2 2 
P + 2 Pq( u + V + w ) = const. 


( 2 . 1 ) 


so that substantial values of v, w imply that p must vary. 


2.3 APPLY ONE-DIMENSIONAL HAVE ANALYSIS 

Assume, reasonably, that the flow near the outer boundary is a small 
perturbation of the free stream. Also assume, less reasonably, that the flow 
in a radial direction resembles a one-dimensional flow. Introduce a radial 

coordinate r and a velocity component u r in that direction. These 

assumptions together imply that 

( 57 + a 0 It )(p + p 0 a 0 u r ) = 0 (2 * 2a) 

( !t ~ a 0 !f )(p " “oVr 5 = °- (2 * 2b) 

Now consider steady flow, and consider two grid points as shown in Figure 
1. Connecting b to i through equation (2.2a) leads to 

(P b ~ Pi> = -P 0 a 0 (u rb " u ri>* (2 * 3a) 

Connecting b to the free stream by (2.2b) leads to 



p 0 > ■ p 0 a 0 <u rb ' U r0 ) ' 


(2.3b) 
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These two equations may be solved for p^, u r b» However, they imply that 

p, - p . u , - u . 

^ (2.4) 

p 0 ’ p b u r0 ‘ “rb 

We expect that in the far field both p and u r decay monotonically to 
their free stream values. In that case both fractions in (2.4) should be 
positive, but here we force them to have opposite signs. Careful examination 
of the output from codes using this method reveals this non-monotone behavior 
at the boundary. 

The observation made in this section originates with Dr. Cedric Lytton, 
of the Royal Aircraft Establishment, Farnborough, United Kingdom. 


3. THE NEW PROPOSALS 

3.1. SPECIFY ENTROPY AT INFLOW AND OUTFLOW 

This is not of course a new proposal. It is perhaps the only widespread 
current practice that is truly unobjectionable. For it to be valid, we merely 
have to draw the outer boundary far enough away that no shockwaves intersect 
it. Since entropy is constant along particle paths in smooth flows, we should 
specify s = sq at points of inflow, and extrapolate s from the interior at 
points of outflow. 


3.2 SPECIFY VORTICITY AT INFLOW AND OUTFLOW 

This is commonly done where vorticity is used to formulate the problem 
(e.g., vorticity-streamfunction treatment of Navier-Stokes equations). The 
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author does not know of its use for solving the Euler equations in conserva- 
tion form. It allows us to specify the tangential velocity components on the 
boundary, equal to their free stream values at inflow or extrapolated from the 
interior at outflow. Probably it is sufficient to make a zeroth order extra- 
polation, but the convection laws of vorticity 



(u»V)m + 


(u)*V )u 


Vp x Vp 

2 

P 


could be used instead. 

In either case, the numerical implementation of this condition should 
probably be through the expression for the integrated vorticity in a cell 


i / m dv = i / u x n ds. (3.1) 

Then in a finite-volume analysis, say, there will be boundary cells with one 
exterior face. All other faces carry fluxes determined by the interior 
scheme, so that enforcing (3.1) will determine the tangential component(s) of 
velocity on the exterior face. Thus we have (n - 1) boundary conditions in 
n dimensions. 


3.3 USE BICHARACTERISTIC ANALYSIS ON THE ACOUSTIC WAVES 

(a) The Two-Dimensional Case 

We seek solutions of the Euler equations at large distances, assuming an 


expansion of the form 
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p = p Q + ep 1 + ... 
P = P 0 + ep L + ••• 


u 


+ eu, 


+ • . . 


v = £Vj + . . . 


Insertion of these expansions into the complete equations gives 


3p ^ 3p ^ 3U| 3vj 

TT + u 0 3x“ + P 0 ( ^T + ^y _) 


= 0 


Su l x 3u l 1 3p l 

Ft - + u 0 F5T + p^ TjT 


3vj avj L apj 

FT + u o ST - + 97" 


= o 


= 0 


9 9 p ^ 2 ^ u i 3 v | 

ft + u o air + p o a o ( ir~ + w } 


= 0 . 


(3.2a) 


(3.2b) 


(3.2c) 


(3. 2d) 


These equations are partly uncoupled, in that p^ does not appear in the 
last three but could be found after these have been solved. From now on we 
exclude the first equation. An equation holding in a characteristic plane can 
be obtained by multiplying (3.2b) and (3.2c) by p^a^cose, p^a^sinO 
respectively and adding them to (3. 2d). The result is (we have dropped the 
first-order suffixes) 
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t|t + (u o + a o cose) + a o sin8 If] 

+ p 0 a 0 cos9[|^ + (u Q + a o sec0) |£] (3.3) 

+ p o a o sine tfr + u 0 + a O cosec0 If] = °* 

It is easy to check that the differential operators inside each bracket all 
act in one plane, which is the property that qualifies (3.3) as a 
characteristic equation. A revealing rearrangement, however, is the following 


8 8 3 

[gY + (u 0 + a Q cos0 ) + a Q sin0 jj] {p + P o a Q (ucos0 + vsin0)} 


2 8 8 

+ Pq3q [ sin0 — - cos0 gy] (usin0 - vcos0). 


(3.4) 


Here the first operator acts along a particular bicharacteristic (TP in Figure 
2) on the sum of pressure plus p Q a 0 ti mes the component of velocity in the 
direction 0. The second operator acts only in space, perpendicularly to 
the direction 0, on the velocity component in its own direction (PQ). 
Writing u r for the velocity in direction SP, and u 0 for the velocity in 
direction PQ, we have 


[ (u 0 ) Q - (u Q ) p ] 

(p + P 0 a Q u r ) p - (p + P 0 a 0 u r ) T + P 0 a Q ^ = 0. 


Even in a steady flow, this does not reduce to equation (2.2a). To employ 
such a relationship as a boundary condition, we need to choose a specific 
value of 0. There are various tempting choices, but we may take a lead 


0 


from the work of Bayliss and Turkel [3J. They showed that under the 


transformation 


n = y, £ = j , t = Ba 0 t + M Q f 


(3.5) 


where M Q = u Q /a 0 and g = (1 - M^) , the system (3.2) implies that p 

obeys a regular wave equation 


P “Pr!-~P = 0 

t t K ££ nn 


(3.6) 


or, equivalently, 


P tt " P RR R P R ~ r 2 P H ° 


(3.7) 


where 


„2 _ 2 2 x 4, . 2 

R =? +n =__+y 

e 


„ . n 3y 

tandi = — = — . 
£ x 


Now at large distances, the last term in (3.7) tends to be small. Foi 
example, if p is given by the well-known separable solution 


p(t , R, 4> ) = e cos(n<|> )J^(kR) 


the orders of successive terms in (3.7) are R ^ R R ^/2^ R 5/2^ 


Therefore, we truncate (3.7) to 
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p — p — — p = 0 

a T h rr r h r 


(3.9) 


with relative error R - ^. We observe that the trial solution 


p = f(r - R)/R 


1/2 


(3.10) 


which represents a decaying outgoing solution, satisfies (3.9) up to the order 
of the neglected term. Also, (3.10) satisfies exactly 


P t + P R 


+ fe " °* 


(3.11) 


Note that this equation holds along an outgoing bicharacteristic of 
(3.9). (This seems a little strange at first, but considering the analogous 
analysis of the one-dimensional case in Appendix A makes it seem natural). In 
fact, (3.11) is selecting, out of all the local bicharacteristics at a point, 
that one which coincides with a global bicharacteristic (see Figure 3). For 
our purposes, that is the most useful bicharacteristic because along it we can 
write an equation derived from global considerations. 

Under the transformation inverse to (3.5), the bicharacteristics of (3.9) 
must become the bicharacteristics of the system (3.2). The distinguished 
bicharacteristic equation (3.11) becomes 


Sp 

Tt 


+ 


BR 


- M q x 




+ y 


(p - p 0 )] = o 


(3.12) 


which is one form of the boundary condition recommended by Bayliss and Turkel 
[3] to suppress incoming radiation. Here also, we take (3.12) to be the 
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equation which will determine pressure on the boundary. The differential 
operator in (3.12) coincides with the bicharacteristic operator in (3.3) if we 
choose 


sin0 




$R - M q x 


(3.13) 


Other useful forms of this result are 


COS0 


(m q x + 8R)y 

sln9 2 x 2 

x + y 

* ■ BM 0 E BxR - V 2 
SR -V ' * 2 + y 2 ‘ 


With this choice of 0 , equation (3.2) can be written 


8 P + 6 U . y ?P) 

at gR - m qX L ax y ay J 


x - pm q r 
+ p 0 a 0 gR - M Q x 



B 3 R a 0 
x - gM Q R 


3u-| 

3x-l 


g 2 y 

P 0 a 0 gR - 


M 0 x 


■3v 

■at 


+ u 


3v 
0 3x 


(gR - M 0 x)a 0 


e 2 y 


ay 


•] = 0. 


(3.14) 


Combining (3.3) and (3.14) leads to 


, ow „ \ 3u _3_ 3u .2 r 3v 3v-i 

(X - BM 0 R) + 8 R »„ 57 + 6 ytjf + u 0 jj] 


+ a o« R - V> §7 ■ 


(3.15) 
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Since (3.12) and (3.14) both hold along an outgoing bicharacteristic, it 
should be a stable numerical procedure to evaluate the spatial derivatives 
from inside the boundary in both cases; hence in (3.15) also. 

Equation (3.15) can be used to update the boundary value of the linear 
combination 

(x - 0 M q R)u + 0 2 yv. (3.16) 


Since our proposal is to use vorticity to update the component of velocity 
tangential to the boundary, we need to ensure that these two conditions are 
independent. In other words, the boundary contour must never lie in the 
direction 


dy = 

dx 


e 2 y 


- fm 0 R 


f (M q , y/x). 


(3.17) 


These prohibited directions are sketched in Figure 5. In the limit Mq+0 

these directions are radial; as Mq+I they are horizontal in the left half- 
plane, in the right half-plane they are tangential to circles that are 
centered on the y-axis and pass through the origin. Clearly, there is little 
temptation to construct any boundary curve that follows these directions. 


(b) The Three-dimensional Case 

No new ideas are involved here, but the procedure is harder to 
visualize. However, we can simply repeat the formalism of the two-dimensional 
case, adjoining to (3.2) the additional equation 
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3w 3w 3p. 

3T" +U 0 3T" + p^3 ^“ !=0 

and taking the three-dimensional divergence in (3. 2d). 

Then a characteristic combination is 

tit - + (u 0 + a 0 cos9) 57 + a 0 sin 9 cos< l > 57 + a Q sine sin<|> |£] 

+ p o a o cos0 \ t £ + (u o + a o sec9) 57] 

(3.18) 

+ P 0 a 0 sin 9 cos <l> [|y + u o 57 + a Q cosec 0 sec<f> |y] 

. „ . r 3w 3w . 3 Wn 

+ p^^sin© sir*J> + u q 97 + a gCOsec0 cosec<|> y^-J . 

It may be checked that all four operators act within one three-dimensional 
space, which is defined by the equation 

cos 0 + sin 0 cos<j> + ^ sin 0 sin<j> = Uq + a^cos© . (3.19) 

Visualized in the space (x/t, y/t, z/t) this is a plane surface (see 

Figure 6 ) which touches a sphere whose radius is ag and whose center is at 

(ug, 0, 0). The line OP is the direction along which pressure is 

differentiated in (3.18). It is bicharacteristic in the sense that it is the 
intersection of planes having parameters (0 ± d 0 ) , (<|> ± d<J>). 

Again, we will select specific values of 0 , <j> by appeal to the far 
field analysis. The transformation (3.5) (with £ = z) produces 


k 
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P - p__ - p ~p =0. 

tt F nn cc 


(3.20) 


At large distances, solutions of this equation have the form 


p = f(i - R)/R 


(3.21) 


where 

R 2 « ? 2 + n 2 + z 2 = * + y 2 + z 2 . (3.22) 

6 

These solutions obey the differential equation 


* P E + E - °- 


(3.23) 


As before, this can be transformed back into an equation along the 
bicharacteristic, allowing pressure to be updated thus 


3p 

Ft 


+ 


e 2 a 0 

3R - M q x 




Z H + P ' P 0 ] = °‘ 


(3.24) 


The differential operator here acts along the bicharacteristic in (3.18) if we 
choose as before 


COS0 


x - gM Q R 
(3R - M q x 


(3.25) 


and also 


tand) = — . 

y 


(3.26) 
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Using these expressions, we convert (3.18) into 


e 2 a 

3 p p u r 3 p 3 p 3 pi 

3t + gR - M q x 3 x + y 3y + Z 3 


x - gM rt R 


3 

S R a r 


- f3u M “0 3u-j 

P 0 a 0 gR - M.x L 3 t x - gM n R 3x J 


g 2 y r 3 v , 3 v , (8R " M 0 x)a 0 3v, 

p 0 a 0 g R - M n x Ut u 0 3x + .2 3yJ 

0 g y J 


(3.27) 


A g 2 z r 3 w 3w (0R “ M 0 x)a 0 3w-, 

P 0 a 0 gR- M n x [<Tt + u 0 lx + ~2 TP 

U p z 


Now we combine the two bicharacteristic equations (3.24) and (3.27) to 


<* - • V > [*T + kI 


0 ~ ~ (SR - M x) a 

, q 2 r 3v 3v 0 3v-i 

** ~2 a 0 5y] 

g y } 

0 a . (SR - M n x) a 

, _2 roW , dW . 0 dWi 


r 0 w i ow, 

+ 8 z t5T + u 0 + 


‘0 


(3.28) 


P - Pr 


This equation allows us to update the velocity component 


(x - gM n R)u + g (yv + zw) . 
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The condition that avoids having this component parallel to the boundary is 
the same as in two dimensions. 


4. DISCUSSION 

The intention announced in the introduction has been carried out. A set 
of boundary conditions, sufficient in number to determine the flow, have been 
obtained from physical considerations. In all cases, it is possible to 
represent the conditions as finite-difference formulae involving the proper 
domains of dependence. Nevertheless, there is a possible hazard associated 
with expressing the inflow boundary condition in terms of vanishing 
vorticity. Effectively this is a derivative condition which fails to 
communicate what the velocity vector at infinity actually is. Although the 
statement is true, it is incomplete. In fact, if the initial data for the 
problem is close to uniform flow, our boundary conditions take the form of 
specifying that no outside influence creates any changes. With a conservative 
scheme, total momentum within the computational domain will change only 
through boundary effects which have been allowed for. What might happen is a 
slow drift away from the desired velocity, which could probably be stabilized 
by prescribing a constant velocity magnitude, say, at one point. 
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APPENDIX A 


The One-dimensional Case 

One -dimensional acoustic flow is governed by the pair of equations 


p t + p 0 a 0 U x = 0 


(A. 1 ) 


u + — p =0 
t P Q x 


(A. 2) 


from which we may deduce the wave equation 


p tt - a 0 p x* ' 0 


(A. 3) 


with its general solution 


p = f(x - a Q t) + g ( x + a Q t). 


(A. 4) 


If there are no incoming waves at large x, then g = 0, and p 
satisfies 

P t + a 0 P x = °* (A,5) 


This corresponds to equations (3.11) or (3.23) in the text. We also have 
the characteristic equations 


4t + a o ?7 )(p + p o a o u) 

(A. 6) 

it " a o §7 )(p “ p o a o u) * 

(A. 7) 


-18- 


At the outer boundary, we discard (A.7), which should carry no information 
inward, and retain (A. 6) which propagates the inner solution outward. 
Combining (A. 6) with (A. 5) yields 

u_ + a_u = 0 (A. 8) 
tux 

so that (A. 5) and (A. 6) are two outgoing characteristic equations, stating 
that p and u remain constant on lines dx = agdt. Thus we see that the 
assumption of no incoming waves enables us to write two outgoing 
characteristic equations, one for p and one for u. This is precisely what 
happens in the main body of the text. 

To test these ideas in one dimension, a simple code was written to solve 
(A. 1 ) , (A. 2) on a grid (0,1,..., N + 1). At points 1 through N, the 
solution was updated by a first-order upwind scheme based on the 
characteristic variables (p ± Pgapu). At 0, N + 1 both p and u were 
updated by first order upwind schemes using data from (0,1) or (N,N + 1). The 
initial data consisted of an internal disturbance superposed on a uniform 
state. The disturbance passed cleanly and stably through the boundary with no 
reflections whatsoever. 
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Figure 5. The forbidden boundary directions 
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the sphere at P 


Figure 6. Geometry of three-dimensional bicharacteristics 
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